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Introduction 

Crystals grown from liquid solutions have important industrial applications. Zeolites, for 
instance, a class of crystalline aluminosilicate materials, form the backbone of the chemical 
process industry worldwide, as they are used as adsorbents and catalysts [Meier, 1986; Barrer, 
1986]. Many of the phenomena associated with crystal growth processes are not well understood 
due to complex microscopic and macroscopic interactions. Microgravity could help elucidate 
these phenomena and allow the control of defect locations, concentration, as well as size of 
crystals. Microgravity in an orbiting spacecraft could help isolate the possible effects of natural 
convection (which affects defect formation) and minimize sedimentation. In addition, crystals 
will stay essentially suspended in the nutrient pool under a diffusion-limited growth condition. 
This is expected to promote larger crystals by allowing a longer residence time in a high- 
concentration nutrient field. Among other factors, the crystal size distribution depends on the 
nucleation rate and crystallization. These two are also related to the “gel” 
polymerization/depolymerization rate. Macroscopic bulk mass and flow transport and specially 
gravity, force the crystals down to the bottom of the reactor, thus forming a sedimentation layer. 
In this layer, the growth rate of the crystals slows down as crystals compete for a limited amount 
of nutrients. The macroscopic transport phenomena under certain conditions can, however, 
enhance the nutrient supply and therefore, accelerate crystal growth. 

Several zeolite experiments have been performed in space with mixed results. The results from 
our laboratory have indicated an enhancement in size of 30 to 70 percent compared to the best 
ground based controls, and a reduction of lattice defects in many of the space grown crystals 
[Sand et al., 1987; Sacco et al., 1993], Such experiments are difficult to interpret, and cannot be 
easily used to derive empirical or other laws since many physical parameters are simultaneously 
involved in the process. At the same time, however, there is increased urgency to develop such 
an understanding in order to more accurately quantify the process. 
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In order to better understand the results obtained from our prior space experiments, and design 
future experiments, a detailed fluid dynamic model simulating the crystal growth mechanism is 
required. This will not only add to the fundamental knowledge on the crystallization of zeolites, 
but also be useful in predicting the limits of size and growth of these important industrial 
materials. Our objective is to develop macro/microscopic theoretical and computational models 
to study the effect of transport phenomena in the growth of crystals grown in solutions. Our 
effort has concentrated so far in the development of separate macroscopic and microscopic 
models. The major highlights of our accomplishments are described bellow. 

Macroscopic Model 

The crystal growth process is modeled as a multiphase (solid and fluid phases) system. The 
macroscopic fluid dynamics model accounts for the continuous solid and fluid phases. 
Throughout the process there is continuous interaction between the two phases. For the time 

being, the crystal growth process is modeled by an exponential law, expressed as r = ae ~ p ' , 
where a and p depend on microscopic fluid and chemical interactions. In the future, this model 
will be combined with a microscopic model based on a Monte-Carlo approach. Following 
[Druzhinin, 1995], the mass and momentum equations are formulated as follows. 


Solid Phase 

The continuity equation is 

| ^+V-(j« v ) = r v (1) 

where s is solid volume fraction, u s is solid velocity vector, r s is rate of crystal growth divided 


by the density of the solid p s . i.e., r 
The momentum equation is 



d(su Xi ) 
dt 


+ V • n v ) = sF, +sg t , 


i = 1,2,3. 


( 2 ) 


where F] is the interaction force between the two phases, and, g, the gravity acceleration 
component each along the coordinate direction, x i . 


Fluid Phase 

The continuity equation is 

+ = (3) 

where u f is the fluid velocity vector and r f the fluid depletion rate divided by the density of the 


fluid, i.e., r s - 



. Notice that r /? =r.p,. The negative sign in r f shows that fluid is depleted 


by supplying nutrients for the crystal growth. 
The momentum equation is 


^ ^ + v ■ ((/ - s)u, u, ) - -(V/p), + (1 - j)s + (v •?),,. = 1,2,3. (4) 
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where 8 is the density ratio p y p . p is the fluid pressure divided by p f . The two phases are 

coupled through the interaction force F In our model, p s and p f are assumed to be constant, r 

is the viscous stress tensor divided by the density of the fluid defined for the incompressible 
Newtonian fluid phase as: 


r tJ =(J-s)v 


du 


du, ') 

+ 


dx , dx. 


i,j = 1,2,3. 


(4a) 


In the above expression v is the kinematic viscosity of the fluid phase. 

Constitutive Relations 

In the case of small heavy particles, the major contribution to the interaction force is from the 
viscous drag force [Maxey and Riley, 1983], For particles with Reynolds number of order of 
unity, the viscous drag force is expressed by the following empirical formula: 

F <=~r( u f, ~ u M Re r) ( 5 ) 


wn d ' 5 

Where r „ = 

" 18v 


, Re v = f — - d, and f(Re p ) = 1 + 0.15 Re 
The following relations are used to describe the growth rate 
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Boundary and Initial Conditions 

On the interior surface of the reactor, a no-slip boundary condition ( u, = 0 ) for the fluid phase 

and a slip boundary condition (u f n = 0) for crystal particle phase is imposed, here n is the unit 
normal vector on the surface. Initially, both phases are assumed to be stationary. The particle 
phase is initialized by a specified solid volume fraction distribution in the reactor. The 
subsequent motion of the interactive phases is induced by gravity. 

Results and Discussions 

Fig. 1(a) shows the solid velocity field at an intermediate time in the settling process of particles. 
As a result of the gravity, the particles settle to the bottom and collect around the center due to 
the action of fluid flow. Fig. 1(b) shows the fluid velocity field whose main feature is the 
recirculation region consisting of two large vortices. Fig. 1(c) represents the fluid pressure 
contour showing that the initial linear static pressure distribution changes mainly in the particle 
region. Fig. 1(d) illustrates the evolution of solid volume fraction along the centerline of the 
container: as time progresses, the particles shift to the bottom being slightly diffused. The 
maximum value decreases with time. Fig. 2. represents the effect of variable gravity levels: the 
larger the gravity level, the faster the particles settle. Fig. 3 shows the particle sedimentation on 
the bottom of the reactor: with larger gravity level, the particles settle faster. 

Microscopic Model 

The model is based on the Direct Simulation Monte Carlo (DSMC) method and is aimed at 
capturing the effects of macroscopic flow on the growth and size distribution of the growing 
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crystals. The system include nutrients, crystals and the mean motion of the fluid in a domain that 
models realistic systems. The grid is based on the mean-free path (T WI ) of collisions between 
nutrients. The length of this grid may be modified to include collisional interactions between 
nutrients and the fluid phase of the solution. Nutrients and crystals move and collide under the 
actin of gravity and fluid-drag forces. 

Nucleation Process 

In a zeolite solution nucleation and crystallization overlap and interact with each other. These 
processes have been modeled in a variety of methods [Thompson and Dryer, 1985]. The 
nucleation will be introduced with a Monte Carlo model that resembles closely the actual 
process. The nucleation sites will be created in random with a probability based on a model that 
describes the nucleation as a function of local thermodynamic. Another alternative is to introduce 
nucleation via a chemically reacting process based on a rate coefficient that depends on 
thermodynamic parameters. Ultimately, this issue will be resolved using experimental data. 

Collisions between Nutrients 

Elastic collisions between nutrients are modeled via the No-Time-Counter (NTC) method [Bird, 
1994]. Collisional pairs are selected in cells and the collision probability in each computational 
cell is calculated based on the NTC method. The model employs the hard sphere (VHS) in 
evaluating collision cross sections. 

Nutrient Impingement and Crystal Growth Process 

The nucleation sites grow to crystals of finite size. Particle methods are ideal to model the 
crystallization process as well. The level of modeling however, is dictated by the computational 
requirements. Our approach is to describe the key physical processes with a model that will be 
refined through successive iterations with data comparisons. For this reason, we initially 
concentrate our efforts in the determination of the size distribution of the crystals. 

The crystals are considered as spherical particles that grow after their attraction and attachment 
of nutrients. Within the framework of the DSMC, an attachment probability P a will be evaluated 
in those cells in which nucleation sites or crystals exist. This probability will be calculated on 
concentrations or surface kinetics based on impingement rates. 

Particle Weighting Scheme 

Given that the crystal particles in our simulations are real particles, different particle weights 
have to be assigned to the crystal and nutrient ‘species’ respectively. The primary concern in 
enacting this scheme is the treatment of collisions between simulation particles that are 
represented by different particle weights. In this event, linear momentum and energy are not 
conserved by the conventional NTC elastic collision model. This is accomplished in our model 
using a conservative multi-weight method. 

Results and Discussion 

Preliminary simulations were performed in order to test the motion and collision algorithms of 
the code. The simulations compared the effects of gravity on settling time of crystals in the 
absence of fluid/particle drag. The domain is initialized with 912 spherical crystal particles in a 
cubic container with size 0.3 m. The spherical crystals have a diameter of 4. 17x1 O' 6 m, and mass 
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of 6.63x1 O' 14 kg. Results are shown in Figures 4 and 5 for 1-g and 0-g conditions respectively at 

0.6 seconds after initialization. At this point, the center of gravity of the crystal particles was 
calculated to be at a height of 0.09 m and 0.073 m respectively for the 0-g and 1-g cases. As 
shown, the imposition of a gravitational force field has immediate effects on particle distribution. 
Future simulations will determine the extent to which this distribution gradient affects the growth 
rate and subsequent final size of the crystals as they compete for nutrients. 
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(d) Evolution of solid volume fraction 
along the centerline. 

Fia. 1 Results at time 0.37 s. for a 2-D reactor, v = 1 0' 6 m e s\ a = 1 .0 m-s' 2 . d = 1 0' 4 m. 
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Fig. 3 Particle sedimentation on the 
bottom forv = 10" 6 m*s , d = 10 m. 
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Fig. 4: Crystal particle distribution after 0.6 
Seconds (0-g) 
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Fig. 5: Crystal particle distribution after 0.6 
seconds (1-g) 
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